---
title: "Output example"
author: "Thomas de Graaff"
date: "October 19, 2015"
output: pdf_document
---

## Introduction

This progress report gives output for all crime of the whole population

## Preamble 

First, we need to read in packages, my own functions and the two datasets

```{r results='hide', message=FALSE, warning=FALSE}
setwd("C:/Users/tgf200/Dropbox/Thomas/papers/Crime")
  library("dplyr")
  library("tidyr")
  library("foreign")
  library("ggplot2")
  library("rootSolve")
  library("AER")
  library("reshape2")
  library("quantreg")
  library("zoo")
  source("./prog/R/iterationBayer.R")
  source("./prog/R/CharacteristicsEq.R")
  source("./prog/R/MakeFig.R")
  source("./prog/R/FindEquilibria.R")

cr <- "crime"
  ####################################################
  # Choose whether estimation for only the youth
  ####################################################  
  youth <- 0
  ####################################################
  # Choose whether only for municipality averages
  ####################################################   
  mun <- 0
  ####################################################
  # Read and manipulate data (still manual selection!)
  ####################################################
  datatemp <- 0
    if (youth) {
    data <- read.csv("./Data/Thomas_data_PC4_crime_youth.csv", header=TRUE, sep = ",")
  } else {
    data <- read.csv("./Data/Thomas_data_PC4_crime.csv", header=TRUE, sep = ",")
  }
  data <- data %>%  # fill in crime type
            mutate(pfield = pfieldcrime, 
                   interaction = pfield * addrdens,
                   alpha = alpha_crime,
                   se = sealpha_crime
                   ) %>%
                   filter(!is.na(pfield))
  dataindividual <- read.dta(paste0("./Data/hat_any",cr,"2006.dta"))
  dataindividual_j <- read.dta(paste0("./Data/hat_any",cr,"2006_jongeren.dta"))
  data <- data %>%
    group_by(gemcode) %>%
    mutate(
      onepermean=weighted.mean(oneperdens, tot_bev, na.rm = TRUE),
      oneparentmean=weighted.mean(oneparentdens, tot_bev, na.rm = TRUE),
      perperhhmean=weighted.mean(perperhh, tot_bev, na.rm = TRUE),
      educationmean=weighted.mean(opleiding, tot_bev, na.rm = TRUE),
      socclassmean=weighted.mean(socklasse, tot_bev, na.rm = TRUE),
      twoearnmean=weighted.mean(k_tweeverd, tot_bev, na.rm = TRUE),
      outmigmean=weighted.mean(v_uit_perc, tot_bev, na.rm = TRUE),
      inmigmean=weighted.mean(v_in_perc, tot_bev, na.rm = TRUE),
      houseownmean=weighted.mean(perchouseown, tot_bev, na.rm = TRUE),	
      polavailmean=weighted.mean(polavail_mean_2005, tot_bev, na.rm = TRUE)
    )
  ### New dataset to be used for quantile regression, keep only the missing values
  data_total <- data %>% filter(is.na(alpha))
  data_total$alpha <- na.fill(data_total$alpha,-3.5)
  data_total$se <- na.fill(data_total$se, 1)  
```

## We then specify the specifications: 

```{r results='hide', message=FALSE, warning=FALSE}
  forminit  <- alpha~addrdens + oneperdens +oneparentdens +
    perperhh + opleiding + socklasse + k_tweeverd + 
    v_uit_perc + v_in_perc + schooldens + perchouseown + shops + polavail_mean_2005+pfield+interaction
  formcrime <- alpha~addrdens + oneperdens + oneparentdens+
    perperhh + opleiding + socklasse + k_tweeverd + 
    v_uit_perc + v_in_perc + schooldens + perchouseown + shops + polavail_mean_2005+pfield+interaction |
    addrdens + oneperdens + oneparentdens +
    perperhh + opleiding + socklasse + k_tweeverd + v_uit_perc + v_in_perc + 
    schooldens + perchouseown + shops + polavail_mean_2005+instrument+instrinter
  forminitmun <- alpha~addrdens + schooldens + shops + onepermean + oneparentmean + perperhhmean + 
    educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean +
    polavailmean+pfield+interaction
  formcrimemun <- alpha ~ addrdens + schooldens + shops + onepermean + oneparentmean + perperhhmean + 
    educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean + polavailmean + pfield+interaction|
    addrdens + schooldens + shops + onepermean + oneparentmean + perperhhmean + 
    educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean + polavailmean + instrument + instrinter
  
  formhelprq1 <- pfield~addrdens + oneperdens +  oneparentdens +
    perperhh + opleiding + socklasse + k_tweeverd + 
    v_uit_perc + v_in_perc + schooldens + perchouseown + shops + polavail_mean_2005 + instrument + instrinter
  formhelprq2 <- interaction~addrdens + oneperdens + oneparentdens +
    perperhh + opleiding + socklasse + k_tweeverd + 
    v_uit_perc + v_in_perc + schooldens + perchouseown + shops + polavail_mean_2005+ instrument + instrinter  
  formrq <- alpha~addrdens + oneperdens +  oneparentdens +
    perperhh + opleiding + socklasse + k_tweeverd + 
    v_uit_perc + v_in_perc + schooldens + perchouseown + shops + polavail_mean_2005+
    pfield+interaction + poly(v1,4) + poly(v2,4)    
  
  formhelprq1mun <- pfield~addrdens + schooldens + shops + onepermean + oneparentmean +
    perperhhmean + educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean +
    polavailmean + instrument + instrinter
  formhelprq2mun <- interaction~addrdens + schooldens + shops + onepermean + oneparentmean +
    perperhhmean + educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean +
    polavailmean + instrument + instrinter  
  formrqmun <- alpha~addrdens + schooldens + shops + onepermean + oneparentmean +
    perperhhmean + educationmean  + socclassmean + twoearnmean + outmigmean + inmigmean  + houseownmean +
    polavailmean+pfield+interaction + poly(v1,4) + poly(v2,4)
  data_total<- select(data_total, pc4, alpha, se, addrdens, oneperdens, oneparentdens, 
                  perperhh, opleiding, 
                   socklasse,k_tweeverd, v_uit_perc, v_in_perc,
                   schooldens, perchouseown, shops, polavail_mean_2005, pfield, interaction, 
                   onepermean, oneparentmean, perperhhmean,  
                   educationmean, socclassmean, twoearnmean, outmigmean, inmigmean, houseownmean,
                   polavailmean)
  data <- select(data, pc4, alpha, se, addrdens, oneperdens, oneparentdens, 
                    perperhh, opleiding, 
                    socklasse,k_tweeverd, v_uit_perc, v_in_perc,
                    schooldens, perchouseown, shops, polavail_mean_2005, pfield, interaction, 
                    onepermean, oneparentmean, perperhhmean,  
                    educationmean, socclassmean, twoearnmean, outmigmean, inmigmean, houseownmean,
                    polavailmean)
  dataindividual$directions.foreign <- factor(dataindividual$foreign)
  dataindividual$foreign <- as.numeric(dataindividual$directions.foreign) - 1
  dataindividual_j$directions.foreign <- factor(dataindividual_j$foreign)
  dataindividual_j$foreign <- as.numeric(dataindividual_j$directions.foreign) - 1
```

## And then the estimation procedure:

```{r}
      output <- iteration2sls(dataindividual, data, data_total, formcrime, forminit, formhelprq1, formhelprq2, formrq, youth=FALSE, initvalue=0)
      summary(output$iv)
```

## With final analysis of the output

We first want to find all equilibria for all iterations

```{r}
matrices <- findequilibria(output)
```

Then we want to find the percentages of 3 equilibra occuring per iteration

```{r}
counteq(matrices$cmat)
```

Then we want to know the number of low equiblibra (smaller than 50%). This also indicates the number of equilibria changing from low to high (larger than 50%)

```{r}
counteqlow(output$instrument)
```

And finally, we want to know whether our found equilibria are close (in this case the difference should be smaller than 2.5% in an absolute sence) to the real crime rates

```{r}
percclose(datatemp$pfield, output$instrument, 0.025)
```
## Figure equilibria

And finally, we end with a figure of the equilibria, which for this case is not very exiting.
```{r}
makefig(output)
```


